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ABSTRACT 

This paper presents a mixed basis approach for Laplace 
eigenvalue problems, which treats the boundary as a perturba- 
^ tion of the free Laplace operator The method separates the 
boundary from the volume via a generic function that can be 



pre-calculated and thereby effectively reduces the complexity 
of the problem to a calculation over the surface. Several eigen- 
values are retrieved simultaneously. The method is applied 
to several 2 dimensional geometries with Neumann boundary 
conditions. 



NOMENCLATURE 

D Diffusion constant. 

A Laplace operator/Laplace matrix. 

A- Eigenvalue of A 



A/7 Free Laplace operator. 
S Surface operator 

<I>(x) Fundamental solution of the Laplace equation. 

0, (jc) Surface charge distribution. 

<I>, (x) Solution to the charge distribution o,-. 

I', Eigenfunction to Af. 



INTRODUCTION 

In physics and chemistry many situations with non-trivial 
boundary conditions give rise to a Laplace eigenvalue prob- 
lem that is not analytically solvable. In such situations nu- 
merical methods, such as finite element methods (FEM) IT], 
fast multi-pole methods (EMM) analytic element methods 
(AEM) |]3], boundary element methods (BEM) and bound- 
ary approximation methods (BAM) ||5| (among many others) 
are used. In this paper an alternative method is presented that 
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uses a perturbation expressed in a mixed basis which effec- 
tively reduces the complexity of the problem to a calculation 
over the surface. The mixed basis approach shares similarities 
with AEM,BEM and BAM, which also are formulated on the 
boundary, but has the advantage of not involving a fully pop- 
ulated matrix over the surface and results in an approximation 
of several eigenvalues of the Laplace spectrum. The method is 
applied on diffusion in 2 dimensions for various domains. 

GENERAL CONSIDERATIONS 

We let / be a real valued twice-differentiable function in 
an «-dimensional (Euclidean) domain. The diffusion equation 
is then defined as 

a,v = DAv. (1) 

The general approach of solving the diffusion equation is to 
express the solution as a linear combination (a Schmidt de- 
composition) of terms that separate into a time-dependent and 
a spatially dependent part v(f,x) = Jli'^ifi{^)8iit)- The spatial 
modes are defined by the (Laplace) eigenproblem 

Afi^l,/, (2) 

where is the separative constant. The collection of all con- 
stants Xj that solve this equation give rise to an eigenvalue 
spectrum of the operator A. This spectrum depends on the 
boundary conditions used and on the shape of the boundary. 
In this paper we use Neumann boundary conditions 

Vfi{dD.)-n = (3) 

where dO. denotes the boundary of the domain and n its nor- 
mal. With Neumann conditions the Laplace operator spectrum 
(also called the Neumann spectrum) consists of non-positive 
eigenvalues. In the absence of the boundary conditions we de- 
note the free Laplace operator by Ap and recall that the eigen- 
functions to Ap are harmonic functions. 



S = Af-A, (4) 

where (as before) Ap denotes the free Laplacian. The norm of 
the operator S is not small, and therefore standard perturbation 
theory does not work. 

In this paper we show that perturbation theory can still be 
useful, but that we need to find a more nontrivial set of basis 
functions. It is clear that, if the obstructions do not completely 
block the diffusion, the large scale structure of the low fre- 
quency modes of A are dominated by eigenmodes of the free 
diffusion operator. More technically, we can say that for wave- 
lengths significantly larger than the size of the obstruction, the 
obstructed diffusion operator recovers diffusive behavior, with 
an altered diffusion constant (i.e. shifted eigenvalues). We 
conclude that the eigenfunctions of the free Laplace opera- 
tor should be part of the basis and that this part captures the 
large scale behavior of the eigenproblem. Furthermore we note 
that these eigenfunctions are known analytically and therefore 
there is no need to actually compute them. 

To find the complementary basis that captures the small 
scale structure we use an adiabatic assumption. This is mo- 
tivated by the fast relaxation of local variations during the 
diffusion process. We treat the local variations close to the 
surface as fast degrees of freedom that can be treated using 
adiabatic elimination, i.e. we assume that the fast degrees of 
freedom equilibrate and are determined by the slow degrees of 
freedom. The local adiabatic equilibrium of A is defined by 
Laplace equation AO = with a nontrivial boundary condition 
defined by the low frequency eigenmodes. In principle each 
low frequency eigenmode of the free Laplacian can be modi- 
fied by adding an appropriate linear combination of solutions 
to Laplace equation so that the resulting function satisfies the 
von Neumann boundary condition. Here however we just add 
different solutions to the Laplace equation to the set of basis 
functions and use perturbation theory to calculate an approxi- 
mate spectrum of the obstructed diffusion operator. The fun- 
damental solution to the Laplace equation solves the problem 



SURFACE EXPANSION 

Our goal is to find an efficient computational tool to ap- 
proximate part of the spectrum defined in Eq. |2] A naive ap- 
proach would be to treat the obstructing surface as a perturba- 
tion of the free Laplace operator, and use standard perturbation 
theory to derive the spectrum of the perturbed operator A. Un- 
fortunately this does not work because the surface obstruction 
cannot be expanded in any small parameter In fact we can 
define a surface operator as 



A4>(x)=8(xo), (5) 

where 5 is a Dirac delta function. In two dimensions the fun- 
damental solution is log|x — xo| and in three dimensions it 
is |x — xo|^'. A general solution of Laplace equation with 
some arbitrary (smooth) boundary condition can be expressed 
in terms of integrals over fundamental solutions of different 
charge or dipole distributions over the surface, e.g. 
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4>,(x) = / <iia,(i) log |x-s| 



(6) 



for a charge distribution o, in two dimensions defined on a sur- 
face (one dimensional boundary). The charge distribution is in 
general determined so that Eq.|6]match the boundary condition 
Here we generate basis functions from a set of smooth 
charge distributions. 

In this paper we are only considering two dimensional ge- 
ometries, so the obstructions have a one dimensional boundary 
that is straight forward to parametrize. The charge distribu- 
tions are defined by low frequency Fourier modes, defined as 
functions of a isometric parametrization on the surface of the 
obstructions. Furthermore, by placing the charges at a small 
distance outside the obstructing boundary and then mirroring 
a charge with opposite sign inside the boundary, we effectively 
create dipole charge distributions that are used to generate the 
basis functions. Parametrization of the obstructing surface in 
the three dimensional case it is more complicated but eigen- 
modes of the surface diffusion operator in Eq.|4]can be used. 

Using the above construction we have defined a mixed 
basis consisting of solutions to the Laplace equation <t>, and 
eigenfunctions to the free Laplace operator ^j. To use per- 
turbation theory the basis must be orthogonal. Since Ap is a 
self-adjoint operator, the free eigenfunctions can easily be con- 
structed to be orthogonal to each other with respect to a scalar 
product over the entire volume. The basis from the Laplace 
equation must however be orthogonalized both with respect to 
themselves and to the free eigenfunctions. Using Eq.|5]we can 
express the scalar product between the free eigenfunction and 
the solutions to Laplace equation as a surface integral 



J^dy^'i(y)a%y). (7) 



Using Eq. |6] the scalar product between two solutions to 
Laplace equation can be expressed as 



/ dx<^iix)^j{x) = / dsids2 0i{si)aj{s2)Q.i\si - S2\), 

JV JSxS 



(8) 



where 



The fact that £1 only depend on the distance between si and S2 
can be realized through a rotation x x' where si and S2 are 
placed on a generic, e.g. horizontal, axis at distance |5 — 
This rotation does not change the volume element dx = dx' . A 
key observation is that £2 is a generic function and can be pre- 
calculated and used in all the scalar products and also for dif- 
ferent geometries. This means that the full orthogonalization 
of the mixed basis can be achieved using only surface integrals. 

Let <t>J denote linear combinations of <I>, and ^ j, such that 
<t>J is orthonormal internally and to the free basis ^ j. Then 
perturbation matrix A, based on the orthogonalized mixed ba- 
sis, has elements on the form /y dx^jls^j = 8,], J dx OjA'Py, 
and / dx <t>jA<I>^. The elements can be expressed as surface 
integrals using similar tricks as for the scalar products: 



ii(|ii-i2|)= / dxlog\si-x\log\s2-x\. (9) 
Jv 



dx {x)A^>j (x) = Jds Oi {s)^>J (s) (10) 
/ dx<i>i{x)A^j{x) ^ dsids20i{si)aj{s2)\og\ SI-S0.1) 

Jv JSxS 



The eigenvalues of the matrix A approximate the spectrum 
of A in Eq. |2]in the interval defined by the eigenvalues of the 
free eigenfunctions. If we use a relatively small number of ba- 
sis functions (typically less than 1000), then the diagonaliza- 
tion of A is very fast. In fact, the most computationally expen- 
sive step is the derivation of the elements in A, and especially 
the double surface integrals appearing in the scalar products 
involving two solutions to Laplace equation. To improve the 
computational costs we could possibly use the fact that the in- 
tegral in Eq. [TT]is identical to the total potential energy from 
two interacting charge distributions that do not interact inter- 
nally. There exists several standard methods for efficient cal- 
culation of potential fields of this type, e.g. multi-pole meth- 
ods 111 and multi-grid methods ||6l. Using these methods it 
should be possible to reduce the cost of calculating the double 
surface integrals to o{ns\og{ns)), where its is the number of 
surface elements. 



RESULTS AND DISCUSSION 

To test the surface expansion, the eigenvalue spectrum of 
the mixed matrix A was calculated and compared with the orig- 
inal spectrum of the Laplace operator A for three different 2 di- 
mensional geometries. The first example is a circle in a square 
domain with periodic boundary conditions. The domain con- 
sists of 160000 elements and A was created with the use of 100 
free eigenvectors and 40 surface eigenvectors. Figure [T] shows 
the resulting spectrum of the first 60 eigenvalues in compari- 
son with an ordinary perturbation calculation where S is treated 
as a perturbation of Ap, and the correct spectrum obtained by 
diagonalizing A. The resultant spectrum of A coincides well 
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with the original spectrum of A and it is notable how poorly 
the naive perturbation captures the spectrum. The second ex- 
ample is a more complicated obstruction, a maze, which has a 
relatively large surface to volume ratio. The structure consists 
of 94249 elements. Figure |2] shows the domain and the spec- 
tra of Af, A and A. The spectrum of A was constructed from 
80 surface eigenvectors and between 100 and 200 free eigen- 
vectors. The figure includes much useful information. The 
first eigenvalues of A coincides well with the eigenvalues of A, 
which implies that the basis of A yields a good estimate of the 
low frequency part of the spectrum. Further up in the spec- 
trum the eigenvalues are not accurately estimated, but with the 
use of more free eigenvectors the basis of A is better capturing 
the spectrum of A. This observation indicate that the commu- 
nication between free eigenvectors, via the surface, is locally 
confined for the low frequency part of the spectrum. The inter- 
nal communication of eigenvectors can be measured by inves- 
tigating the eigenvectors of A and proposes an error estimate 
by measuring the weight of the associated eigenvectors of A, 
where few elements of an eigenvector to A implies a good es- 
timate of the eigenvalue. This has however not been studied 
in detail. Finally we present a randomly generated domain 
of 40000 elements. In this domain a large number of surface 
eigenvectors was needed to achieve good results, probably due 
to the large number of sharp corners present. Figure |3] shows 
the domain and the spectrum of Af, A and A, where 250 sur- 
face eigenvectors and 400 free eigenvectors where used. A 
resonant influence was found in the situations where the wave 
number of the eigenvectors coincided with the geometry stud- 
ied. For such eigenvectors, the corresponding eigenvalues of A 
converge slower to the actual eigenvalues with increasing num- 
ber of free eigenvectors. This can be seen in figure|4]where the 
relative error of the eigenvalues of A are plotted against their 
index for the circular obstruction domain. The peaks corre- 
spond to wave numbers that coincide with the radius of the 
circle. 



SUMMARY AND CONCLUSION 

We have constructed a mixed basis that can be used in 
perturbation calculations of the spectrum of the Laplace oper- 
ator in complicated geometries. This reduces the problem to a 
generic function over the volume, which can be pre-calculated, 
and surface integrals. Relatively few vectors are needed to 
span the perturbation matrix which yields for a computation- 
ally interesting approach of solving problems involving the 
spectrum of the Laplace operator, e.g. diffusion. Existing 
techniques such as FEM, AEM, BEM and BAM gives good 
estimates of the Laplace spectrum but generally involve many 
elements. AEM, BEM and BAM are closely related to our 
method, as the calculations are also formulated on the sur- 
face, but for AEM and BEM; require a full diagonahzation 




i 



Figure 1 . ORDINARY PERTURBATION (DASHED), SURFACE PER- 
TURBATION WITH 100 FREE EIGENVECTORS AND 40 SURFACE 
EIGENVECTORS (SOLID) AND LAPLACE OPERATOR SPECTRUM 
(DOTTED). INSERTED PICTURE (ABOVE) SHOWS STRUCTURE 




Figure 2. FREE SPECTRUM (DASHED), A SURFACE PERTUR- 
BATION WITH 200,300 FREE EIGENVECTORS AND 80 SURFACE 
EIGENVECTORS (SOLID) AND LAPLACE OPERATOR SPECTRUM 
(DOTTED). INSERTED PICTURE (ABOVE) SHOWS STRUCTURE 
AND (BELOW) ZOOM OF THE FIRST 50 EIGENVALUES. 

over the surface, and for BAM; an iterative scheme, which is 
needed for each eigenvalue. The mixed basis approach gives 
the opportunity of reducing the size of the resultant matrix to 
derive an approximation to a part of the spectrum. We have 
demonstrated that the resulting spectrum is a good estimate 
and that relatively few elements are needed for good results. 
We have not yet compared the performance of the new method 
to existing techniques. A possible application of our method 
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could be to initialize the solution scheme for an iterative eigen- 
value/eigenfunction solver, for example the method by Li et 
al Q that relies on an initial guess on the eigenvalues. 
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Figure 3. FREE SPECTRUM (DASHED), SURFACE PERTURBA- 
TION WITH 400 FREE EIGENVECTORS AND 250 SURFACE EIGEN- 
VECTORS (SOLID) AND LAPLACE OPERATOR SPECTRUM (DOT- 
TED). INSERTED PICTURE (ABOVE) SHOWS STRUCTURE AND 
(BELOW) ZOOM OF THE FIRST 50 EIGENVALUES. 
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Figure 4. RELATIVE ERROR OF EIGENVALUES OF A OF MIXED 
BASIS APPROACH, COMPARED TO REAL SPECTRUM OF A. A 
WAS CONSTRUCTED WITH 40 SURFACE EIGENVECTORS AND 
THE NUMBER OF FREE EIGENVECTORS WAS VARIED FROM 100 
TO 200 TO 300. INSERTED PICTURE SHOWS STRUCTURE. 
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